Quasi-stationary states and the range of pair interactions 
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Abstract 

"Quasi-stationary" states are approximately time-independent out of equilibrium states which have 
been observed in a variety of systems of particles interacting by long-range interactions. We in- 
vestigate here the conditions of their occurrence for a generic pair interaction V{r- — >■ oo) ~ l/r~' 
with 7 > 0, in d > 1 dimensions. We generalize analytic calculations known for gravity in d = 3 to 
determine the scaling parametric dependences of their relaxation rates due to two body collisions, 
and report extensive numerical simulations testing their validity. Our results lead to the conclusion 
that, for 7 < d — 1, the existence of quasi-stationary states is ensured by the large distance behavior 
of the interaction alone, while for 7 > d — 1 it is conditioned on the short distance properties of the 
interaction, requiring the presence of a sufficiently large soft core in the interaction potential. 

PACS numbers: 05.70.-y, 05.45.-a, 04.40.-b 



In recent years there has been renewed interest in the 
statistical physics of long-range interactions (for a review, 
see e.g. [3), a subject which has been treated otherwise 
mostly in the astrophysical literature for the specific case 
of gravity. The defining property of such interactions 
is the non-additivity of the potential energy of a uni- 
form system, which corresponds to the non-integrability 
at large distances of the associated pair interaction, i.e., 
a pair interaction V{r — >■ 00) ~ with 7 < d in d 

space dimensions. The equilibrium thermodynamic anal- 
ysis of these systems is very different to the canonical 
one for short-ranged interactions (with 7 > d), leading 
notably to inhomogeneous equilibria as well as other un- 
usual properties — e.g. non-equivalence of the statistical 
ensembles, negative specific heat in the microcanonical 
ensemble. Studies of simple toy models have shown that, 
like for gravity in d = 3, these equilibria (when defined) 
are reached only on time scales which are extremely long 
compared to those characteristic of the mean field dy- 
namics. On the latter time scales one observes typically 
the formation, through "violent relaxation" , of so-called 
"quasi-stationary" states (QSS), interpreted theoretically 
as stable stationary states of the Vlasov equation (which 
describes the kinetics in the mean field limit). In this 
letter we consider whether the occurrence of such QSS 
driven by mean-field dynamics can be considered as a 
behavior arising generically when there are long-range 
interactions in play. Using both simple analytical results 
and numerical simulations, we argue for the conclusion 
that it is only for 7 < d — 1, i.e. when the pair force is 
absolutely integrable at large separations, that QSS can 
be expected to occur independently of the short distance 
properties of the interaction. For 7 > d — 1, on the other 



hand, their occurrence will be conditioned strongly also 
on short distance properties, and thus cannot be con- 
sidered to be a result simply of the long-range nature 
of the interaction. Our analysis shows the relevance of 
a classification of the range of interactions according to 
the convergence properties of forces rather than potential 
energies which has been formalized in [^j. 

We proceed by first generalizing a calculation origi- 
nally given by Chandrasekhar for Newtonian gravity to 
a system of N particles interacting by a pair potential 
V{r) = (where g is a coupling constant). This calcu- 
lation, which numerical studies indicate is accurate both 
parametrically and quantitatively for gravity (see e.g. 0- 
01), will give us an estimate of the relaxation rate due 
to two body collisions (i.e. the inverse of the time-scale 
on which a typical particle's velocity is randomized by 
such interactions). Denoting by Tmf the characteristic 
time for the formation of a QSS (i.e. of the mean-field 
dynamics), the criterion for the existence of QSS we will 
then study is 



^2 Tmf — !■ when N 00 . 



(1) 



where the limit N 00 corresponds to the mean-field or 
Vlasov limit . Indeed if this condition is not satisfied, 
it implies that there is no mean-field regime in which QSS 
may form. 

Following the treatment for the case of gravity (see e.g. 
3], section 1.2.1) we consider a test particle of velocity v 
crossing a system in a QSS, assumed spherical and of ra- 
dius R and approximated as homogeneous. We estimate 
first the rate of relaxation due to soft two-body collisions 
by calculating Au^ , the mean square velocity change of a 
particle per crossing (i.e. in a time of order Tmf) due to 
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such collisions. It is straightforward to show that 



dx 



X 



(2) 



where bmin is the minimal impact parameter at which 
the scattering is soft (i.e. the deflection angle is small), 
defined by 



\9\ 



(3) 



and bmax is the maximal impact parameter for two body 
collisions. In these formulae, and in what follows below, 
we use the symbol ~ to indicate that the numerical fac- 
tors in all expressions have been dropped, leaving only 
the parametric dependences which are relevant to our 
considerations here. In the case of gravity in c? = 3 the 
choice of bmax has been a source of debate, with numeri- 
cal simulations indicating that bmax R accounts better 
for results than the more evident choice bmax the 
mean inter- particle separation (see e.g. Q). We consider 
in what follows both possibilities, and will see that our 
central results are not in fact sensitive to which is correct. 
We have also implicitly assumed o? > 1 and 7 > 0. 
We now write 



r2 T'm/ — r^o/t Tm/ + ^hard Tmf , 



(4) 



where the first contribution is that considered above, and 
the second is the remaining one from hard scatterings, 
i.e., collisions with impact factors b < bmin- Taking 
now that Tmf ~ it is straightforward to deduce from 
Eq. ([2]) that, for sufficiently large N, 
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if 7 < (rf- l)/2 
if 7 > (rf- l)/2 



(5) 

if bmin/bmax ^ 1 for large N. To infer these scalings we 
need only (as in the corresponding derivation for the case 
of gravity ^]) use the fact that the QSS is, by definition, 
a virialized state, i.e., we take 
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where U, the total potential energy of the QSS, and K, 
its total kinetic energy, have a fixed ratio because of viri- 
alization. This scaling with N corresponds to that in the 
usual mean field or Vlasov limit, in which U and K both 
scale in the same way with N. 

Using again the scaling Eq. ([SJ, the definition Eq. ^ 
gives 



(7) 



Note firstly that this implies bmin/bmax — >■ as iV —> oo 
for any 7 > if bmax a-nd for any < 7 < c? if 

bmax ~ £ ~ i?7V~^/'*, so that Eq. (O is indeed valid in 



TABLE I: Summary of two body collision rates (without core) 



< 7 < 


soft collisions at ~ bmax dominate 


< 7 < d- 1 


collisions at ~ 6mm dominate 

^ soft Tmf ^ N T yrhardTmf 


7 > d - 1 


Tso/t Tm f and Fhard Tm f divergent in 



these cases. Using now again Eq. ([7|) in Eq. ([5]) we obtain 
the scaling 



^ soft ^mf 



+ if^<(d-i)/2 

^-(d-l-7)/7 if^ > (d- l)/2 



(8) 



where J = if bmax ^ R, and S = {—2j + d — l)/d 
if bmax ^ RN-^I'^. It follows that, for 7 > d - 1, the 
contribution of soft two body scatterings alone diverges 
at large iV, so that the criterion ([Ij cannot be satisfied 
in this case for the "candidate" QSS. For any 7 < d — 1, 
on the other hand, the contribution T soft Tmf vanishes as 
N 00. It is simple to show, in this case, that Thard Tmf 
also goes to zero when N ^ 00, and thus that the condi- 
tion ([1]) for the existence of QSS may be satisfied. To do 
so it is sufficient to consider that this contribution can be 
bounded below by that from an "exactly hard" core with 
radius e = bmin, i-e., U(r) =00 for r < bmin- Estimating 
the collision rate on such a core as The ^ nav where n 
is the mean density and a ^ e'^~^ (the cross section), we 
obtain 



he Tmf 



N 



^-(d-l-7)/7 



(9) 



when we take e = bmin, with the latter scaling as in 

Eq. It follows that T hard Tmf < The Tmf ^> aS 

N oo for 7 < d—1. Further it follows from the inferred 
scaling of TheTmf that, for 7 < d — 1, the total rate T2 
will scale as calculated for Tsoft in Eq. In other 

words, an exact calculation including Thard should give, 
at most, a r2 larger than Tsoft by a numerical factor. 

A corollary of these results, which are summarised 
in Table „ is that, for a QSS to exist in the case that 
7 > d—1, the pair potential must include a suffi- 
ciently large soft core. Indeed to remove the divergence 
of 

Tso/t Tmf in this case, we must introduce a smooth- 
ing of the potential at a scale e which vanishes more 
slowly than bmin in Eq. ([7]). In this case Tsoft Tmf is 
given by the second expression in Eq. ([5]) but with bmin 
replaced by e. Keeping e/R constant, for example, gives 
Fso/t Tmf ^ as iV — )■ cxD for any 7. If the core 

is "exactly soft", i.e., V{r) = g/e"' for r < e we have 
Thard — and the satisfaction of the condition Eq. ([T]) 
follows. If the core is hard, as envisaged above, it is clear 
that the same is not true. Indeed it is simple to check 
using Eq. ^ that it is not possible to choose e in order 
to satisfy both The Tmf — > and Tsoft Tmf simulta- 
neously as cxD for 7 > d — 1. 

These results lead then to the primary conjecture of 
this article: for pair potentials with V{r — >■ 00) ~ ^/r'^, 
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QSS can always exist if there is a sufficiently large soft 
core, but only for 7 < c? — 1 can they exist when such a 
core is not present (i.e. when its size e 0). The validity 
of this conclusion rests evidently on the assumption that 
the dominant correction to the mean-field dynamics is, 
just as for gravity in d = 3, two-body collisionality. More 
specifically we also require the parametric dependences 
of the inferred relaxation rates, which have been derived 
using various simplifying approximations (notably that 
of homogeneity both in configuration and velocity space). 
We now present results of numerical simulations (in d = 
3) which test their validity. We focus here on the crucial 
result above: the parametric dependence of the two-body 
scattering rate due to soft scatterings in Eq. ([S]), for the 
range 7 > (d — l)/2. 

We perform molecular dynamics simulations us- 
ing a version of the publicly available gravity code 
GADGET2^8i]. We have modified the force routine in 
the tree-PM version of the code to treat a generic power- 
law pair potential with a core. As in the original code we 
use a soft repulsive core, with compact support: for r < e 
V{r) decreases continuously to a minimum at r « e/2 and 
then increases back to a local maximum V{r = 0) = 0. 
In what follows the values of e quoted correspond to the 
separation at which the force is still attractive but has 
dropped to approximately 30% of its value in absence 
of smoothing. We consider here the attractive case (i.e. 
g < 0). The simulations are checked using simple con- 
vergence tests on the numerical parameters, and their 
accuracy is monitored using energy conservation. For 
the time-steps used here it is typically of order 0.1% over 
the whole run, orders of magnitude smaller than the typ- 
ical variation of the kinetic or potential energy over the 
same time. As initial conditions we take the N particles 
on the sites of a simple cubic lattice of side Lq and as- 
cribe random velocities uniformly distributed in an inter- 
val [—A, A] in each direction (i.e. "waterbag" type initial 
conditions in phase space). The parameter A is chosen 
so that initial virial ratio is unity, i.e., 2K/\U\ = 7. We 
choose this initial condition because it would be expected 
to be close to a QSS to which (collisionless) relaxation 
should occur "gently" . The system is enclosed in a cubic 
box of side L w 2Lo (and centered on the same point as 
the initial cube of particles). Energy conserving "soft" 
reflecting boundary conditions are used in the dynamics, 
i.e., at each time step particles which have moved outside 
the box have the appropriate components of their veloc- 
ity inverted. The results we report here required runs 
lasting as long as two weeks on up to sixteen processors. 

Shown in Fig. [T] is the evolution of the potential energy 
(7 as a function of time, for a pair potential with 7 = 5/4 
and e/L = 0.01, for the different values of N indicated. 

We have defined t„j/ = y , which, given that 

Lq ~ R, is equivalent parametrically to the deflnition 
used above, assuming the scaling in Eq. ([6]). The macro- 
scopic behavior monitored in this plot is clearly very con- 
sistent with what has been anticipated, in line with the 
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FIG. 1: Temporal evolution of the total potential energy U 
divided by its initial value U{0), for 7 — 1.25 and a soft core 
e/L — 0.01, for the different indicated. 



typical behavior observed in self-gravitating systems and 
other systems with long-range interactions studied in the 
literature: there is a first phase of "violent" (collision- 
less) relaxation towards an approximate equilibrium, the 
QSS, which then evolves itself in a second phase on a 
time-scale which clearly depends on N. The first phase, 
on the other hand, should be A^-independent: as N in- 
creases we see that the different curves are increasingly 
well superimposed at early times. 

Shown in Fig. [2] are, for the two cases 7 = 5/4 and 
7 = 3/2, our measurements of the relaxation rate 

^ relax , 

as a function of A^ (upper panel) at a chosen fixed e, and 
as a function e (lower panel) at fixed chosen TV. The 
estimate of Treiax is obtained simply from the slope of 
the potential energy plotted as a function of time, in the 
region in each case where this is well fit by a linear behav- 
ior, i.e., we take Treiax = d{\n U) / dt at t 0. Each point 
corresponds to one numerical simulation. Note that for 
these determinations we consider thus only the evolution 
away from, but still close to, the QSS. Further results on 
the longer time evolution of these systems, and in par- 
ticular the compatibility of the fully relaxed states with 
those predicted analytically for this case in (and re- 
lated numerical studies in [10]) will be given elsewhere. 

The upper panel of Fig. [2] includes a line showing the 
scaling proportional to 1/A^ predicted by Eq. ([SJ at fixed 
bmin = £■ The agreement is clearly very good. Further 
it is simple to verify that the results are quantitatively 
very coherent with the prediction: taking i? w Lq/2 w 
L/4, Eq. ([5]) fit the normalizations of the plot with a 
pref actor of order unity in both cases. While the degree 
of this concordance — despite the many approximations 
which inevitably limit the accuracy we can expect, and 
the fact that we have dropped all numerical factors in our 
derivation — is clearly just fortuitous, this quantitative 
coherence of the results confirms their solidity. 
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FIG. 2: Estimated relaxation rate rmiax as a function of A'^ 
at fixed e/L = 0.01 (upper panel), and as a function of e at 
fixed A'^ — (lower panel). 



The lower panel of Fig. [5] shows likewise excellent 
agreement with the predictions above. On it are shown 
lines corresponding to the behavior of Eq. ([5]) at fixed N, 
when we replace bmin by e. As discussed above, this scal- 
ing is predicted to be valid in the regime bmin < e < bmax ■ 
Below bmin we expect the rates to reach an asymptotic 
e- independent value of order those estimated in Eq. ([5]). 
The behaviors in the plot are very coherent with these 
predictions, for values of bmin which are in good agree- 
ment with Eq. d?]), taking again Rk, L/A. While we have 
not predicted the value of bmax, the downward deviation 



(corresponding to a reduction in scattering rate) from the 
fit at larger e occurs at a value very consistent in each 
case with the measured mean inter-particle distance. 

We note that this last plot explains why we consider 
only 7 up to 7 = 3/2, and indeed why we have not tried 
to verify more directly the scalings in Eq. using sim- 
ulations with e — ^ 0. The reason is that, in order to mea- 
sure the relaxation rates, we need to access the regime 
^reiaxTmf 1, i-G., we need to have a reasonable separa- 
tion between the times scale of the collisionless dynamics 
(and formation of QSS!) and the relaxation time scale. At 
N = 10^ we see that 7 = 3/2 is already at this limit for 
the smallest e, and the error bars on these points refiect 
the greater difficulty we have in making the measurement 
in these cases. The only remedy is to increase N, which, 
however, is prohibitively expensive numerically, in partic- 
ular at small e where the proper integration of the (few) 
hard collisions included requires significant decrease in 
the time stepping. 

Finally a few remarks on the relation of these results 
to some of the extensive recent literature on QSS (see [H 
for references). The determination of the N dependence 
of QSS lifetimes has been much emphasized, both as a 
target for phenomenological studies of toy models, and 
for theoretical studies of the problem. Our results show 
that such lifetimes can be expected to depend, in general, 
not just on N, but also on the parameters characterizing 
the short distance properties of the potential. While for 
7 < d — 1 a limit e = may be defined [and gives the 
scaling of Eq. ([5])], for 7 > d — 1 this is not possible and 
the scaling of the relaxation rate will depend necessarily 
on how e scales with N. It would be interesting to ex- 
tend our numerical simulations to explore the robustness 
of QSS notably to effects which may come into play in 
more physically realistic settings: as shown by a recent 
study of a toy model, the introduction of stochastic- 
ity in the dynamics may also destroy QSS. We emphasize 
that our results here apply only to the particular (albeit 
broad) class of models considered, and a priori not, e.g.., 
to long-range spin models in which there is no equivalent 
of two body collisions. It remains an interesting open 
question to determine in a broader such context the con- 
ditions for the existence of QSS on the spatial dependence 
of the interaction. 

We acknowledge many useful discussions with F. 
Sicard. This work was partly supported by the ANR 
09-JCJC-009401 INTERLOP project. 
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